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A new approach to cosmological perturbation theory has been recently introduced by Bartelmann 
et al., relying on nonequilibrium statistical theory of classical particles, and treating the gravita¬ 
tional interaction perturbatively. They compute analytic expressions for the nonlinear matter power 
spectrum, to first order in the interaction, and at one-loop order in the linear power spectrum. The 
resulting power spectrum is well behaved even at large wavenumbers and seems in good agreement 
with results from numerical simulations. In this paper, we rederive their results concisely with a 
different approach, starting from the implicit integral solution to particle trajectories. We derive the 
matter power spectrum to first order in the interaction, but to arbitrary order in the linear power 
spectrum, from which the one-loop result follows. We also show that standard linear perturbation 
theory can only be recovered at infinite order in the gravitational interaction. At finite order in 
the interaction, we find that the linear power spectrum is systematically and significantly underes¬ 
timated. A comprehensive study of the convergence of the theory with the order of the interaction 
for nonlinear scales will be the subject of future work. 


I. INTRODUCTION 

The statistical properties of cold dark matter (CDM) 
in the nonlinear regime make for a technically challeng¬ 
ing problem, and their study has been the bread and 
butter of several generations of cosmologists. Various 
approaches have been pursued to obtain analytic predic¬ 
tions for the matter power spectrum and higher-order 
statistics beyond the well-established linear perturbation 
theory results [T]. The reader may refer for example to 
Ref. [2] for a detailed review and to Ref. [3] for a list or 
more recent works. 

The most standard analytic approaches to large-scale 
structure rely on an expansion in variables measuring 
departures from perfect homogeneity. The perturbation 
variables are the density and velocity field in the case of 
Eulerian perturbation theory (see e.g. Ref. [J]), or the 
displacement held in the case of Lagrangian perturba¬ 
tion theory (LPT, see e.g. Ref. [5]). In a given gauge, 
all perturbation schemes agree on a unique linear den¬ 
sity held, which matches observations on large scales, on 
which the departures from homogeneity are indeed small. 
On quasilinear scales, perturbative approaches improve 
upon linear theory, but all methods eventually fail dra¬ 
matically at small enough scales, where inhomogeneities 
cannot be treated perturbatively [H]. 

In a series of recent papers im, Bartelmann et 
al. (hereafter BFB) have proposed a new approach, based 
on non-equilibrium statistical theory of classical particles 
and using tools similar to the path-integral approach of 
quantum held theory. Their approach introduces a new 
expansion parameter, in addition to (and independently 
of) the usual perturbation variable^ the order in the 
gravitational interaction. To set our lexicon right away. 


^ In BFB the expansion with respect to the usual perturbation 
variables is carried out when evaluating the initial phase-space 
probability distribution. The latter is computed perturbatively 


throughout this paper we will use “linear”, “quadratic”, 
“one-loop” or “non-linear” to refer to the dependence on 
standard perturbation variables. When discussing the 
interaction expansion, we shall say explicitly “to n-th or¬ 
der in the interaction”. BFB derive analytic formulae 
for the nonlinear power spectrum to hrst order in the in¬ 
teraction, which are well behaved at large wavenumbers 
and seem to match results from numerical simulations 
deep into the non-linear regime, way beyond the reach of 
standard perturbation schemes. 

In this paper we rederive all of BFB’s results with 
a completely different approach, starting from the im¬ 
plicit integral solution for particle trajectories. We ob¬ 
tain explicit solutions for the trajectories or equivalently 
for the CDM phase-space density, perturbatively in the 
gravitational interaction. The density field and its power 
spectrum are then explicitly computed to first order in 
the gravitational interaction. Our solutions are on the 
other hand non-perturbative in the standard perturba¬ 
tion variables (or in the initial departure from inhomo¬ 
geneity), therefore going beyond BFB’s results, which are 
restricted to one-loop order in the linear power spectrum. 
From these nonperturbative solutions, we show how BFB 
obtain their approximate solutions. We argue that they 
rely on an improper treatment of the “damping factor”, 
which leads to a spurious enhancement of power at small 
scales, mimicking the effect of nonlinear growth. 

We also make contact with standard linear perturba¬ 
tion theory. We show that it can only be recovered at 
infinite order in the gravitational interaction, and with a 
rather slow rate of convergence as a function of the order 
in the interaction. This warrants a careful study of the 
convergence of the theory in the nonlinear regime, where 
explicit calculations so far (including in the present work) 


in the small initial correlations between particles’ positions and 
momenta, which are proportional to the initial matter power 
spectrum. See Appendix A of Ref. [3- 
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have been limited to first order in the interaction. 

This paper is organized as follows. In Section |ll] we 
recall standard results about clustering in the Zeldovich 
approximation, from which the results of Ref. [5] are im¬ 
mediately derived. In Section III we define the new ex¬ 
pansion scheme and rederive the main results of Ref. [5] 
starting from particle trajectories. Section |TV] makes the 
connection with standard linear perturbation theory. We 
conclude in Section Ivl 


II. DARK MATTER STATISTICS IN THE 
ZELDOVICH APPROXIMATION 



A. Summary of known results 

In this section we re-derive the results of Ref. [8] , where 
the statistics of CDM in the Zeldovich approximation are 
computed approximately. 

The well-known Zeldovich approximation is the first 
order of LPT [5] . The statistics of CDM in this approx¬ 
imation have been extensively studied, see for example 
Refs. [TUIITT] . 

In this approximation particle trajectories are given by 

= q + ip{q,ri), (I) 

where x is the comoving coordinate, q is the conformal 
time, q is the unperturbed Lagrangian coordinate and 
is the displacement field. The latter is a curl-free gaus- 
sian random field. It is related to the linear density field 
i5l through V • '4>{x,q) = —SL{x,q), or equivalently, in 
Fouriei[^ space 

%l){k,q) =-i^5-L{k,q). (2) 

On subhorizon scales, the linear density field grows with 
a scale-independent growth factor: 


??) = D{q)6i,{k, qi), (3) 

where qi is some initial time. This implies that the dis¬ 
placement field remains parallel to itself and that tra¬ 
jectories are straight lines. In particular, the peculiar 
velocity field is 


dX'Zel 

dq 




(4) 


The Zeldovich approximation is only accurate to lin¬ 
ear order in the density field, but one can still compute 
the clustering statistics that would result if particles fol¬ 
lowed these trajectories. The overdensity field S{x) is 


^ Our Fourier transform convention is 
f{k) = fd^x e*'“ “’/(a:), f{x) = fd^k/{2nf e-‘'“ “’/(fc). 


FIG. 1. Variance of the displacement field per axis, per log¬ 
arithmic fe-interval, at redshift zero: = ^kPi,{k)/{2^^). 

The linear power spectrum was obtained with Camb |13| for 
cosmological parameters consistent with Planck results |14 |. 


then obtained by summing over the streams leading to 
the position x: 


lp5ie\{x,q) = J d^q ( 5 d (x 

where 5 d is the Dirac distribution, 
transform, we get [TU] 


-Xzei{q,q)), (5) 
Taking the Fourier 


Szeiik) = J d^q - I 


( 6 ) 


The power spectrum [defined by {6{k)6*{k')) = 

(27r)^P(fc)5D(fe)] follows immediately [I^ : 


Pze\{k) = e 


-k^ 


d^q e*''-'* - ij , (7) 


where Cij{q) = {ipi{0)tj}j{q)) is the correlation tensor of 
the displacement field and is given by 

In equation Q, Pl is the power spectrum of the linear 
density field, and is the variance of the displacement 
field (per axis). It is obtained from 


2 1 f d^k PL(fc) 

3 7 (27r)3 fc2 


(9) 


At zero separation, we have Cij{0) — cr'^dij. 

We show the characteristic variance of the displace¬ 
ment field per axis, as a function of scale, in Fig. 
We see that the Zeldovich displacement field fluctuates 
mostly on scales k <1 Mpc“^. 

The main advantage of the Zeldovich approximation 
over linear theory is its better correlation with the fully 
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nonlinear density field m This arises because the Zel- 
dovich approximation accounts for the advection of struc¬ 
ture by large-scale flows, which cannot be accounted for 
by standard linear perturbation theory m- 

It is however well known that the nonlinear power spec¬ 
trum in the Zeldovich approximation is actually lower 
than the linear prediction [12) . and a fortiori lower than 
the nonlinear power-spectrum. Indeed, in this approx¬ 
imation, structure is exponentially suppressed on scales 
smaller than the characteristic displacement, tr,/, « 9 Mpc 
at redshift 0 HD, i.e. for wavenumbers fc ^ ~ O.I 

Mpc“^. This is a consequence of the unphysical as¬ 
sumption that particles move forever on straight lines, 
even after streams cross. Actual particle trajectories are 
strongly deflected at stream crossing, and CDM eventu¬ 
ally forms bound structures. To account for this requires 
a framework to treat nonlinear evolution, using for ex¬ 
ample higher-order analytic methods or numerical simu¬ 
lations. 

We now Taylor-expand equation Q to second order in 
Pl (recalling that tr^ is of the same order as Pl), and 
arrive at [T^ 


Pzeiik) « PLik) - k^alPLik) + Pi.ioop(fc), ( 10 ) 


where 


^l-loop(^) 


/k-k'Y fk- {k-k')\ 

^ (vW; V (fc-feT ) 


( 11 ) 


As expected (and by construction!), the linear term is 
just the linear power spectrum. The quadratic terms 
arise from the nonlinear mapping between the density 
field and the displacement field. The term —k'^a^Ph{k) 
cancels out the contributions of the long-wavelength 
modes k' k oi the loop integral to the power spec¬ 
trum P{k). Physically, it embodies the fact that as long- 
wavelength modes of the displacement field cannot have 
an effect on the small-scale power spectrum [16) . 

Since the Zeldovich approximation is only the first or¬ 
der of LPT, the resulting power spectrum is only mean¬ 
ingful and accurate to first order in the linear power spec¬ 
trum, and should only be trusted for k < « 0.1 

Mpc“^. Obtaining the correct term of order P^ in the 
CDM power spectrum would require going to third order 
in perturbation theory. 


B. Discussion of BFB’s result 


we shall see) recover the known Zeldovich power spec¬ 
trum, Eq. Q. 

The statistical treatment of BFB requires them to 
compute the joint probability distribution of the posi¬ 
tions and momenta of N particles at some early ini¬ 
tial time. They partially expand this probability dis¬ 
tribution, to second order in the position-position and 
position-momentum correlations, but keeping the full ex¬ 
ponential dependence on the correlation matrix of mo¬ 
menta (Sections IV E and IV F of Ref. [7]). As a re¬ 
sult, they arrive at a partially expanded power spectrum 
for the Zeldovich approximation (Eqs. (44) and (48) of 
Ref. [5]). With our notatioij^ it takes the form 

Pzel,BFB(A:) = e~^ [PL(fc) + Pl-loop(fc)] , (12) 

where the 1-loop term is given by our equation 0- At 
this point, BFB correctly remark (Section V B of Ref. [5]) 
that since the term in brackets only accounts for corre¬ 
lations of the initial phase-space distribution to second 
order, it is necessary to also approximate the exponential 
damping prefactor at the same order. However, instead 
of writing a consistent Taylor expansion, which would 
lead to our Eq. ( |Io| ) , they approximate the damping term 
e~^ by (l-|-fc^cr^)“^, and apply it only to the one-loop 
term, so their final result is 

^Zel,BFB,final (^) = Phik) + ^ ^ ' (1^) 

1 -k 


Their justification of this choice is that the exponential 
suppression in the Zeldovich approximation is unphys¬ 
ical, and should therefore not be applied to the linear 
term. The choice of replacing the negative exponential 
by a fraction is not discussed. 

For fc > 1 Mpc“^, the loop integral is converged when 
including only fc' <C fc and \k' — k\ ^ k contributions (see 
Fig. 0, and is approximately 

^’l-loop(fc > 1 Mpc"^) « fc^cr^PL(fc)- (14) 


This implies that equation Eq. (13) has the asymptotic 
value 


-Pzel,BFB.final(fc ^ 1 MpC ^ 2 X Pi,{k), 


( 15 ) 


which is precisely the behavior seen in Fig. 1 of Ref. [5]. 

Whereas the resulting power spectrum appears to be 
closer to a nonlinear power spectrum than what one 
would obtain from the Zeldovich approximation (it is en¬ 
hanced instead of exponentially damped at small scales), 
arriving at Equation (13) results from a series of choices 
that are neither mathematically consistent nor physically 
justified. 


In Ref. [5] , BFB derive the expected clustering of par¬ 
ticles treated with the “Zeldovich propagator”, i.e. en¬ 
forcing that their trajectories are those given by the Zel¬ 
dovich approximation. Without additional physical in¬ 
gredients to deal with the deflection of trajectories at 
stream crossing, their treatment should (but does not, as 


® The notation correspondence is crJ/3 = H-ri = D(r])/D{r]i). 
Ref. [2 finds a quadratic contribution proportional to instead 
of (1 -b Ti)^. This difference arises from the choice of initial 
conditions: a gaussian density field for BFB, and a gaussian 
displacement field in our case. It is unimportant at late times. 
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C. Bispectrum 

We now compute the bispectrum B{ki,k- 2 ,k^) [de¬ 
fined as usual by {5{ki)5{k2)5{k^)) = (27r)^(5D(fei -t fc 2 -I- 
k 3 )B{ki,k 2 ,k 3 )] in the Zeldovich approximation. We 
recall that the Zeldovich approximation is only the first 
order of LPT and therefore is not of high enough order 
in standard perturbation variables to compute the bispec¬ 
trum. One should indeed go to second order of pertur¬ 
bation theory to compute the correct bispectrum and it 
should not be a surprise that the two differ. With that 
caveat in mind, we first expand the density field Q to 
second order in the displacement: 

Szei{k) « S^ik) - J d^q 

1 f d^k' k ■ k' k ■ {k - k') 
’^ 2 ] {k'Y {k-k'Y 

x5i^{k')6Uk-k'), (16) 

where in the second line we have used Eq. ([^ to express 
ifi in terms of the linear density field. The lowest-order 
bispectrum then follows in a straightforward manner: 

Bzei(fci, fe 2 , fca) = F{ki,k 2 )PL{ki)Pi^{k 2 ) perm., 

(17) 

where 


Denoting the scale factor by a, the conformal time by q 
and the comoving position by a;, the equations of motion 
of a point particle are given by 


dx 

dq 

du 

dq 


u 

7 

a 


—aV(/). 


(19) 

( 20 ) 


The Newtonian gravitational potential (j) satisfies the 
Poisson equation: 


A(/) = 47rGpa^(5 = WqO (21) 


where <5 = p/p — 1 is the overdensity field, which has to 
be determined from the trajectories themselves, and we 
have defined, for short. 


Wg = 47rGp(z = 0) = 


( 22 ) 


The system (19)-(20) has the following implicit integral 
solution: 


xiji) = x{q^) -I- (s - s^)u{q^) 

- f {s-s')a'V(j){x{q'),q')dq', 
dv 

where p* is some initial time and 


(23) 


F{ki,k2) 


(fci -f k2) ■ ki (fci -I- k2) ■ k2 

P P 


. kik2 
kik2 



(fci • fe2)^ 


, ( 18 ) 


is identical to equation (79) of Ref. [5], and different, 
as expected, from the correct result obtained in second 
order (Eulerian or Lagrangian) perturbation theory [2]. 

Note that here we have expanded the density field first, 
for the ease of computation. One can also easily compute 
the bispectrum first, from the unperturbed expression 
for 6{k). When doing so one recovers the exponential 
damping prefactor of Ref. [H] . 

So far we have rederived part of the results of BFB, in 
the illustrative case of the Zeldovich approximation. We 
now turn to the most interesting and novel result, where 
the gravitational interaction is treated perturbatively. 



is the “superconformal time”. 

We assume that up to the initial time p*, particle tra¬ 
jectories are well described by the Zeldovich approxima¬ 
tion, so our initial conditions are: 

x{q*) = q-\-ip{q,q*), (25) 

) 

u{q^) = a., s iPiq^q*) = K'ip{q,q*), (26) 

D{q*) 

with V •'0(q, p*) = —Shiq, p*). For compactness we shall 
define 


a(p) = 1-I-A*(s — s*). (27) 

In what follows we shall write ■0(g,p*) = 'ifiq), the de¬ 
pendence on p* being implicit. 


III. PERTURBATIVE TREATMENT OF THE 

GRAVITATIONAL INTERACTION B. The gravitational interaction as a perturbation 


A. Trajectories of point particles in an expanding 
universe 

The starting point of our derivation are the equations 
of motion of a non-relativistic point particle in an ex¬ 
panding Universe. We shall give a more formal (and 
exactly equivalent) description with the Vlasov-Poisson 
system in Section p^II F| 


A standard technique in field theory is to treat interac¬ 
tions perturbatively. The new conceptual contribution of 
BFB is to use this technique for gravitational interactions 
in cosmology. In this context this approach is somewhat 
artificial as there is no explicit small parameter associ¬ 
ated with the interaction - in fact, as we will show in the 
next section, even the simple linear fluid equations re¬ 
sulting in the linear growth factor are formally of infinite 
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order in the interaction. Nevertheless, it should in prin¬ 
ciple converge to the correct answer as higher and higher 
orders in the interaction are included. It may moreover 
provide insights into the nonlinear regime. 

In order to keep track of the order in the gravitational 
interaction, we insert a bookkeeping perturbation param¬ 
eter e in front of the potential. We shall set e to unity at 
the end of the calculation. 

Our starting point is therefore the following implicit 
solution for particle trajectories: 

x{ri) = q + a{r])ip{q) 

— e / {s — s')a'^(j){x{r]'),q')dr]'. (28) 

•t?). 

We shall expand the trajectories and density field to first 
order in e: 

x{t]) = x^°\q) + e x^^\q) + 0{e^), (29) 

5(a;) = (5(°)(a;) + e ^(^^(a;) + ©(e^). (30) 

The gravitational potential itself, being proportional to <5 
through the Poisson equation, is also a series in e. Since 
it is multiplied by e in Eq. ( [^ , we only need to use 
the zero-th order [sourced by when computing 
trajectories to first order in e. 

We therefore have, explicitly, 

x^°^q) = q + a{q)-ip{q) (31) 

for the free trajectory, and 

= — f (s — s')a''V(j)^^\x^^\q'),T]')dq' (32) 

Jrj, 

for the first-order correction. 


The Fourier transform of V(5 d(®— a^o) is just — 
and we thus get 

5(i)(fe) = i J d\k- a;(i)(q) (37) 

From the Poisson equation, we get 
aV4>^°\x,q) =iujQ J 


rl^y h' , . 


(27r)3 {k')'^ 


(38) 


so that Eq. (32) becomes 


a;(i)(q,) = J J 




Inserting this into Eq. (37), we arrive at 

d^k' k ■ k' 


r'l r h ■ h 

5^'^\k,q)=ujl J dq{s-s') J j^^j^I{k,k',q,q), 

(40) 


where 


I{k,k',q,q')= f 

^(41) 

Defining, for short, K = a{q)k and K' = a{q')k', we 
may rewrite this as 

I(fc,fc',?7,77') = f d^q 5^°\k',q') 

(42) 

As a technical aside, note that (a;) has zero mean, as 
can be seen by interchanging statistical averaging with 
averaging over the whole space in Eq. 




C. Density field to first order in the interaction 

The density field is obtained by summing over all 
streams: 

1-I-(5(a;) = J d^q S-d[x — x^^\q) — e x^^\q)]. (33) 

Expanding to first order in e gives 

1 + A) (a;) = J d^q 5D[a: - a:^°^(9)]. (34) 

and 

6^^\x) = — J d^q x^^'^{q) ■'VSb[x — x^^\q)]. (35) 

Since the free trajectory = q -|- a t/; is, up to the 
coefficient a{q), identical to the Zeldovich trajectory, we 
immediately obtain S^^\k) from equation (1^: 


D. Power spectrum to first order in the interaction 

We are now ready to compute the power spectrum. To 
first order in e, we have 


P{k) = P^°°'>{k) + 2 e P(°i)(fc). 


(43) 


Here is the power spectrum of which again is 

almost identical to the one in the Zeldovich approxima¬ 
tion with the additional coefficient a^\ 

p(oo)(fc) = Jd^qy'^ '^ ^k'ya^CiPq) (44) 

The first-order term comes from the cross-power of 
and It requires computing the following ensemble 

average: 

(I(fe,fc',ry,AA)*(fc",77)) 

= / d^q d^q' d^f^'y(^k~k'yq+rk'-q'-ik"-q" 


S^°'>{k,q) = / d^q e*'“ '» - I 


(36) 


X (e 


J{K~K')-lp{q)+iK' ■tf){q')-iK''-ipiq") 


(45) 
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where K" = a{r])k'' . Note that we got rid of the con¬ 
stant part of which is possible since has zero 

mean. We also got rid of the constant part of 
since the Dirac function Sj^ik') would vanish when mul¬ 
tiplied by k! in the final integral. 

For any linear combination A of gaussian random vari¬ 
ables with zero mearj^ {e^) = e^^^ 1. Working out the 
variance of the last term in brackets, we arrive at 


E. Expansion in Pl 

We now Taylor-expand our result in Pi^{rj^) in order 
to recover the result of Ref. The expansion of 
follows directly from our results for the Zeldovich approx¬ 
imation, equation ( 10 ) with an additional multiplying 
Cij and cr^. In particular, the linear part is 


(...) = exp 


- [{K' - Kf + [K'f + {K"f] 




L°\k) = a'^{r])Pi^{k,r]^). 


(51) 


-{K^-K'^)K'^C,,{q' -q) 

T^ffj / 


+ -K'^)K"^C,j{q" -q) 


+ K'^K"^C,,{q" -q') 


This corresponds to equation (77) of Ref. [3]. 

We now expand The integrand of equation (48) 


(46) 


Changing integration variables to q^qi = q' — Q and 
^2 = q" — q, we see that the integral over q leads to 
(27r)^(5D(fc" ~ fc)i as expected from statistical homogene¬ 
ity. We then get 


{I{k,k',q,q')S^°^*{k",q)) = 
Pxs{o){k,k',q,q'){2TTf6-D(k" - k), 


(47) 


with 


gets expanded into three linear terms and six quadratic 
terms. Only the terms that mix q^ and ^2 survive inte¬ 
gration, which leaves one linear term and four quadratic 
terms. The linear term is 

Pi5m,i.{k,k' ,q,q') 

= aa'{2T:f5j^{k'-k)Py J Sq (q) 

= aaPi,{k, ?7*)(27r)^dD(fe^ — k), (52) 

where we recall that Cij is the correlation tensor of the 
displacement field at 77 *, hence the resulting Pl (j)* ) • The 
linear part of pO^l is therefore 


rv 

PlSi 0 ){k,k\q,q')=e-"^^° Pl°^\k,q)=U}^aiq) dq'[s - s')a{q')Pi^{k,q,). (53) 


exp 


{K'^ - K^)K'^C,j{q^) + - K'^)K^CM 


+ K^K'^C,,{q^-q,) 


where we have defined, as in Ref. [9], 


Qd = 2tT; 




K' -K ■ K' 


(48) 


(49) 


From this we finally obtain the contribution to the power 
spectrum at first order in the interaction 


P'^°^\k,q) =uj^ J dq\s-s')J 


d^k' k k' 


(27r)3 [k')'^ 

X Pj^g(o){k,k',q,q'). (50) 


Equations (43), (44), (48) and (50), along with the inter¬ 


mediate definitions, constitute one of our main results. 
They give the power spectrum to first order in the inter¬ 
action, but to infinite order in the linear power spectrum. 
We have indeed only expanded in e and not in Pl. This 
formal result therefore goes beyond that of Ref. [9] , which 
we shall now derive. 


This is equivalent to equations (78)- (790 of Ref. 0. 

The four quadratic terms of the integrand of Pxsw) 

_( 2 A) _( 2 i^) 

correspond to the terms Z to Z of Ref. [5], given 
in their equations (63) to ( 66 ). For example, the term 
involving the product of Cij{q 2 ) and Cij{q 2 — Qi) is two 
times 

(PT* - 

X j d^qi dPq 2 e*'' ''‘^~''^'’^^Cij{q 2 )Cmn{q 2 - 9i) 

= (PT* - K'")K^K"^K''^ 

X J d^q[ d\2 e-^’^' ‘‘'^+d>^'-'^y‘i-aj{q2)Cmn{q[) 

= (PT* - K'")K^K"^K‘ 


/*^ r^l T^m {k'" - P)(k'^ - y) 


|fc'-fe|4 (P)4 

X Px{k - k',qt)Ph{k',q*) 

K K' K {k- k'){K - K') ■ {k - k') 

{k'Y |fe'-fe|4 

X Px{k - k',qt)Px{k',q*). (54) 


Short proof: if A is a linear combination of gaussian variables, 
A is itself a gaussian variable. Having zero mean, its moments 
are (A^P+i) = 0, (A^p) = (2p)!/(2Pp!) x (A^)^. We then have 

(e^> = Ep(A 2 D/( 2 p)! = Ep((^'>/ 2 )VP! = 


^ The correspondence of notations is gqp{r,T') = A*(s — s'), 1 + 
gqp{T,r^) = a{r]). 

^ Equation (79) of Ref. [9] has a typographic error: the second 
argument of the two gqp should be r* instead of 0 
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_(2A) 

This is the term Z of Ref. [9]. Similarly, the term 

involving and Cij{q 2 — qi) gives and the 

_(2C') 

term involving Cij{qi) and Cij{q 2 ) gives Z . Finally, 
the square of Cij{q 2 — qi) leads to 


where 

Sf{x,u,ri) = f{x,u,ri) - Suiu) (58) 

is the perturbation from a perfectly cold homogeneous 
matter. The initial conditions are 


X J _ q^)Cmn(q2 “ Qi) 

= K^K'^K^K''"{2TrfdB{k' - k) 

J d\ e-*'“ 'JQ, (q)C„„(q) 

r d'^k" 

= {2TTf 5]^{k'- k) J j^^P^ik-k",q^)PL{k'',q^) 

K ■ k" K' -k” K-{k- k'')K' ■ (fe - k") 

(F)4 |fe-fc"|4 ■ ^ ‘ 


Up to a missing factor of (Stt)^, this is the term Z 
of Ref. [5]. Inserting this into equation (50) leads to 
equations (67) and (69) of Ref. [5]. 

Here again, as in the case of the Zeldovich approxi¬ 
mation, Ref. [S] approximate the term e^'^n/z i^y _|_ 
(5 d/ 2)~^, which they only apply to the quadratic term. 
As we pointed out previously, this approximation is not 
correct neither formally nor physically, as the term pro¬ 
portional to —QdPl is missing. In addition, the asymp¬ 
totic behavior at large k of the resulting power spectrum 
can be shown to be proportional to Ph{k). Again, this is 
due to the choice of the approximation for the damping 
prefactor, and the extending of this approximation much 
beyond its regime of validity. 

In order to test the accuracy of the power spec¬ 
trum to first order in the interaction against numeri¬ 
cal simulations, one would need to numerically evaluate 
Eqs. (48) and (50). This involves computing a triple 
three-dimensional integral, which is challenging, and is 
deferred to future work. 


F. Going beyond the first order in the interaction: 
pertnrbed Vlasov-Poisson system 

We can extend our calculation to higher orders in the 
interaction. To do so, we may formulate the problem in 
terms of the Vlasov-Poisson system. In the relevant limit 
of an infinite number of particles, this is exactly equiv¬ 
alent to solving particle trajectories, but takes a more 
compact form for a formal description. 

The Vlasov-Poisson system is comprised of the colli¬ 
sionless Boltzmann equation for f{x,u,q), the phase- 
space distribution (divided by p{z = 0)), and the Poisson 
equation: 

'll 

dnf + - ■ da:f - e aVcj)- d^f ^0, (56) 

a 

A(j) = I (fu6f{x,u,r]), (57) 


f{x,u,r]^)= / (fqSB[x-q-ip{q)]Sr,[u-Xtip{q)]. 

(59) 

Expanding / and (j) in orders of the interaction, we have 
the following explicit recursion scheme: 

= 0, (60) 

a 

=u;^a-^ J dPu (a;, m, r;),(61) 

n — 1 

= a y (62) 

The inhomogeneous partial differential equation 
u 

d^F -\ -• dxF = S(x, M, rj) (63) 

a 

is most easily solved in Fourier space: 

F{k, M, q) = M, 77*) 

-b r '>'^■'‘8{k,u,q'). (64) 

• 71 ). 

The initial condition for the distribution function is, in 
Fourier space, 

/(fc, u, 77*) = J d^q - XMq)], (65) 

so that we have 

f^°'>{k,u,q) = j d\ - A*i/:(q)],(66) 

/(") (fc,«, 77) = -iujI d77'e‘(^-^>'^ J ^ 

n-l . / 

y J("-1-P)(fc',^') . a„/(p)(fc - fc',M,77'). (67) 

p =0 ^ ' 

The first-order calculation would recover the results we 
have obtained directly from particle trajectories. We see 
that going to higher orders in the interaction leads to a 
higher and higher number of terms, as well as a larger 
and larger number of nested integrals, if one is to keep 
the full nonlinear dependence on the initial conditions. 
This quickly becomes intractable for any practical com¬ 
putation. 

A similar iterative solution to the Vlasov-Poisson sys¬ 
tem was suggested in Ref. HZl. There, the zero-th or¬ 
der solution is the phase-space density in the Zeldovich 
approximation, which accurately describes large scales, 
contrary to our zero-th order, non-interacting solution 
(as we shall see in the next section). The hierarchical 
solution of Ref. m is therefore potentially more accu¬ 
rate than an expansion in the gravitational interaction, 
but suffers from the same computational challenges when 
going to higher orders. 









IV. CONTACT WITH STANDARD LINEAR 
PERTURBATION THEORY 


The approach we have followed can allow us in princi¬ 
ple to go to higher orders in the interaction, even though 
explicit computations would become more and more in¬ 
volved. It is best to understand beforehand what hope 
one may have of getting accurate results. This is best 
done by considering the power spectrum on large scales, 
where standard linear perturbation theory is very accu¬ 
rate. 

Since standard perturbation theory does not rely on 
an expansion in the gravitational interaction, it is, for¬ 
mally, of infinite order in the interaction. In this section 
we shall study the convergence of the perturbative inter¬ 
action treatment to standard linear perturbation theory 
on large scales. 

Since we are here concerned with linear theory (specif¬ 
ically, linear in (5 l(^*) or in the initial displacement field 
■ 0 ), we may treat the cold dark matter as an ideal fluid. 
The linearized continuity and momentum equations are 


dr] 


= -0l, 


d 

dr] 


{aOif) 


—a e A(p 


-e Wo(5l, 


( 68 ) 

(69) 


where we have inserted the bookkeeping parameter e in 
front of the gravitational potential and have used the 
Poisson equation in the second equality. This equation 
has an implicit integral solution: 

dhir]) = (5l(??*) + a*dL(?7*)(s - s*) 

+ e ujq f (s - s')Si,{r]')dri'. (70) 

•I 17 , 


If we treat the gravitational interaction perturbatively, 
i.e. consider e as a perturbation parameter, we first write 


(71) 

71 


The subsequent orders in the interaction can then be ob¬ 
tained from the following recursion relation: 

= <^l(??*) + a*^L(?7*)(s - s*) 

= a(r])6L(v*), (72) 

6^\r])=u}l f (s - s')d‘i^~^\r]')dr]', (73) 

dri. 


where a was introduced in Eq. (27). The full solution is 
therefore 


dhiv) 


r 

5h(v ) ^ ^ y (s - s')a(r]')dr]' 

+ ^ (wg)^ / (s - ^)dr( f (s' - s")a(r]")dr]" + ... 

'd n* ^ r}^ 


= T 


>7) 
^0 


V- 

a(r]) 


(74) 


where T is the time-ordering operator. 

We now recall that e = 1 is just a bookkeeping pa¬ 
rameter. Solving the system of fluid equations (|^-(69) 
with e = 1 is equivalent to computing the full solution 
(74) and leads to the usual linear growth factor D(r]): 


d-L(v) = -^^4(?7*)- 
D(v*) 


(75) 


On the other hand, truncating the series (741 at any finite 


order in e does not recover the correct linear growth fac¬ 
tor. Again, we emphasize that this is because standard 
linear perturbation theory accounts for the interaction to 
infinite order, even if it is, on the other hand, only lin¬ 
ear in the initial conditions Si^(r]s,). The two expansion 
schemes (in (5 l(? 7*) and e) are independent and this result 
is perfectly consistent. 

We can now easily recover our previously obtained lin¬ 
ear power spectrum at first order in the interaction: 

PL(v)^PAir]) + 2eP^°^\7]) 


0 ^( 77 )-f 2e Woa(ry) / di]'(s - s')a(r]') Pl(? 7 *), (76) 


which corresponds to equations (51) and (53) obtained 
in the previous section. 

Now how accurate an approximation do different or¬ 
ders in the interaction expansion provide to the correct 
linear growth factor? We first notice that for all n, 
Si^) > 0, as can be seen from the recursion relation ( |73| ). 
Therefore any truncation of the series (7^ will systemat¬ 
ically underestimate the correct growth factor. We show 
in Figure the linear growth factor summed up to order 
n in the interaction, for n = 0,1,2,3 and 00 , the latter 
case corresponding to the standard linear growth factor. 
Our starting time corresponds to redshift z* = 99. We 
see that the series converges very slowly to the correct 
answer at redshift zero, and still grossly underestimates 
the correct linear growth factor even at third order in the 
interaction. 

To compensate for this missing power. Ref. m only 
use the interaction expansion starting from a rather low 
redshiflQ The interaction expansion is indeed more ac¬ 
curate closer to the starting redshift. The choice of the 
starting redshift is however arbitrarily 

In closing of this section, we point out that BFB do 
recover the standard linear growth factor at large scales 


^ Typically around redshift 10 or so m- 

® In Ref. [^, this redshift 2 :* is presumably determined by requir¬ 
ing that the evolution with the Zeldovich approximation from 
some early initial time Zf to 2 *, followed by the perturbative in¬ 
teraction approach from 2 :* to redshift 0, result in the correct 
linear growth factor for large, linear scales. This can never be 
the case since, as we have shown, the perturbative interaction 
approach systematically underestimates the linear growth fac¬ 
tor. It is therefore unclear how one should choose the starting 
time. 
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FIG. 2. Linear density field up to order n in the interaction: 
for n = 0,1, 2, 3 and oo from bottom to top. The normaliza¬ 
tion is such that the linear growth factor (corresponding to 
n < oo) is unity at redshift zero. The starting conformal time 
rjt corresponds to redshift z, =99. Cosmological parameters 
consistent with Planck results are used. 


when they replace the hamiltonian propagator, correctly 
describing particles trajectories, by the “Zeldovich prop¬ 
agator” (equation (65) of Ref. [7]). This amounts to using 
the Zeldovich approximation. They do not derive the Zel¬ 
dovich propagator nor the linear growth factor from first 
principles with the perturbative interaction approach. In 
this section we have explicitly shown how to do so, by go¬ 
ing to infinite order in the interaction. 


V. CONCLUSIONS 

Recently Bartelmann et al. have suggested a new ap¬ 
proach for computing statistics of CDM in the nonlinear 
regime. Their approach is inspired by field theory and 
relies on a formulation similar to the path integral. The 
main physical ingredient of this approach is to treat the 
gravitational interaction perturbatively, as one would do 
for other fundamental interactions in the weakly interact¬ 
ing regime. This approach is perfectly valid from a formal 


point of view and potentially very promising, calling for 
further examination. 

In this paper we have tackled two points, one regarding 
the technique and the second the essence of the approach. 

First, we have concisely rederived all the results of BFB 
with a different formalism, relying on the integral solu¬ 
tion of particle trajectories. We have furthermore derived 
the matter power spectrum to first order in the interac¬ 
tion but to infinite order in the linear power spectrum, 
going beyond the one-loop results of BFB. In addition, 
we have laid out the iterative equations one would have 
to solve to go to higher order in the interaction, using 
the Vlasov-Poisson system. 

Second, we have made contact with standard linear 
perturbation theory and examined the validity of treat¬ 
ing the gravitational interaction perturbatively. We have 
shown that standard linear perturbation theory accounts 
for the gravitational interaction at infinite order (even 
if it accounts for initial conditions at linear order). We 
have also shown that the convergence of a perturbative 
interaction approach is very slow on large, linear scales. 

It is yet unclear whether such an approach would fare 
better in the nonlinear regime. The properties of the 
density field on highly nonlinear scales indeed depend 
on the details of bound or marginally bound structures, 
which may be difficult to reproduce if treating the grav¬ 
itational interaction perturbatively. Perhaps the theory 
can produce accurate results if one applies it repeatedly 
over short enough intervals of time, though one then has 
to account for non-gaussian initial conditions. A rigorous 
study of these issues calls for a detailed analysis of the 
nonlinear regime using the interaction expansion, which 
we defer to future work. 
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